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ABSTRACT 



The advent of small, portable super microcomputers has enabled the rapid analysis of 
mechanical signals. In addition to the microcomputer's reduction in cost, their capabilities 
have expanded to allow real time analysis of the frequency characteristics of vibrating 
structures. Data acquisition modules allow the addition of powerful components to a 
microcomputer. Modules consisting of the necessary Analog to Digital converters, Digital 
to Analog converters, and the hardware necessary to support sixteen channels of input data 
can easily fit inside desk top microcomputers. Fast 32 bit processors, small high volume 
hard disk drives and large RAM capacities complete the package, enabling a microcomputer 
to perform complex dynamic signal analysis. One application of the super microcomputer 
is in the field of underwater explosion shock testing. The combination of the acquisition 
and analysis capabilities in the same device eliminates the requirement to record the data to 
an intermediate device prior to analysis. Thus the data can be analyzed within minutes of 
the completion of an experiment, giving immediate results at the testing site instead of an 
analysis facility. Another application of the data acquisition equipped microcomputer is in 
the analysis of the frequency response characteristics of structures. Programming the 
microcomputer to acquire the data from a vibrating structure in order to determine the 
dynamic signature is easily done. The Hilbert transform, identified previously as a means 
of detecting non-linearities in the frequency response, can also be included in the data 
analysis portion of the same program. Thus a complete package of signal analysis 
routines, including frequency response, power spectra, and system linearity, can be 
programmed into the microcomputer. Discussed in this study is the method utilized to 
validate a data acquisition equipped microcomputer for such applications. 



THESIS DISCLAIMER 



The reader is advised that the computer codes presented in this research have not 
been tested under all possible operating conditions or the various operating systems for the 
MASSCOMP 5400 system. Effort has been made to run the programs on the most current 
operating system available and to the best of the researchers ability in the time available, to 
ensure that the programs are free of computational and logical errors. Any use of these 
programs without user verification is at the risk of the user. 
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I. INTRODUCTION 



A. BACKGROUND 

The super microcomputer is a tool that enables the engineer or scientist to rapidly 
digest large quantities of information. Recent advances in computer technology have 
brought about a startling increase in the computing power of small "desktop" systems. The 
speed of the processor has significantly increased while the cost have dropped dramatically. 
Coupled with the ability to store large quantities of data in internal memory, the speed of 
computers have taken quantum leaps forward. Five years ago random access memory 
banks of 600 kilobytes (600K) were common for a desktop scientific computing system. 
Today systems with 10 Megabytes are readily available. Processor word size has increased 
from 8 and 16 bit word length to 32 bit in the same period. Data storage and retrieval 
systems have also become smaller and less expensive. Thus the ability to process 
information rapidly and store it has made micro and mini computer systems a valuable aid 
to the engineer. 

One of the uses of the super microcomputer is the area of data acquisition. 
Experiments can be rapidly set up, the data can be collected and analyzed in near real time, 
and the results can be displayed in fraction of the time than was previously required. 
Evidence of the increase in computing power is evident in the increased ability of Analog 
To Digital conversion. Analog signals, frequently voltages from a set of remote sensing 
devices, can be automatically read and recorded. In 1973 the sampling rate of 20,000 
samples per second was considered good for a microcomputing system [Ref. 1]. In 
1986 the super microcomputer is capable of sampling 1,000,000 samples per 
second [Ref. 2]. The only limitation is the throughput capability of the data bus. This is a 
50 fold increase in sampling rate and coupled with simultaneous sampling cards and 
multichannel A/D converters, super microcomputers are now able to analyze signals in 
wide band vibration experiments and shock analysis applications. 

In the field of underwater explosion testing the engineer has been forced to record 
shock data to an intermediate device prior to analysis. This device, usually a wide band 
tape recorder, is a costly addition to the test equipment. If there was a small, portable 
device which could record the data and perform the analysis at the test site the engineer 
would have immediate results available. Not only would this save time, but it would 
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eliminate the need for the intermediate recording device. The properly equipped 
microcomputer would allow the engineer to gather the information, display the results, and 
record them for further analysis. This process would take seconds, where as the present 
alternative would take much longer. 

In the study of the mechanical signature of structures, there has been a recent emphasis 
toward use of the Hilbert transform in order to determine if the frequency response of a 
structure is behaving in a non-linear manner. There are no dynamic signal analyzers 
available that include the capability of performing the Hilbert transform of the frequency 
response. In order to perform this operation the data would have to be transferred to a 
computer for further analysis. This is a costly and time consuming process. If a 
microcomputer with a data acquisition package could perform the dynamic signal analysis 
of mechanical signatures from a structure, then it could also perform the Hilbert transform 
of the frequency response of the structure. Again, this would result in a significant 
reduction in the time required for data analysis and include processes that are not available 
on the present generation of dynamic signal analyzers. 

B. PURPOSE 

The purpose of this study is to: 

1 . Program the properly equipped microcomputer to perform data acquisition, reduction, 
and analysis of underwater explosion shock data. 

2. Perform dynamic signal analysis on mechanical signatures from structures, using 
both impact and steady state excitations. 

3. Include in the dynamic signal analysis portion of the data analysis routines the Hilbert 
transform. This would place a complete signal analysis package in one device. 
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H. THEORY 



A. DIGITAL SIGNAL PROCESSING FUNDAMENTALS 

The Fourier transform is a method of calculating the frequency spectrum of a time 

signal. Fourier stated that any periodic time signal of fundamental frequency could be 

*0 

transformed into the summation of a series of sinusoidal signals of varying amplitude and 
frequency. The Fourier series of a periodic time signal of period T 0 is then: 



oo 




n=l 



[Eqn. 1] 



The Fourier transform of a signal is a method of computing the frequency spectrum 
of a non-periodic infinite time signal, provided it satisfies Dirichlet's Conditions. It is 
defined as: 



X(f ) = T{x(t)} = J x(t) e -J 27lft dt . [Eqn. 2] 

-OO 

The finite time signal is infinite in the frequency domain, the infinite time signal has a 
finite bandwidth of frequency components. The fast Fourier transform (FFT) is a recent 
mathematical technique allowing the rapid computation of the frequency components of a 
time signal. Programed into a microcomputer with the ability to convert a time signal into a 
set of digital samples, the FFT is a powerful tool that allows fast computations of a systems 
transfer function. 

The process of digitizing an analog signal is to multiply it by a pulse train at a set 
sampling frequency. Given an analog signal of bandwidth F max , the Nyquist sampling 
frequency states that the sampling rate of the analog signal must be greater than or equal to 
2 F max in order to avoid aliasing errors, thus: 

F s > 2 F max . [Eqn. 3] 
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Multiplication in the time domain is equivalent to convolution in the frequency 
domain. This simply takes the frequency components of the signal and repeats them in the 
frequency domain at intervals of the sampling frequency. If the signals bandwidth exceeds 
the sampling frequency then the higher frequency components overlap and sum to 
introduce errors in the frequency spectrum of the time signal. This is known as aliasing. 
Figures 1 though 3 illustrate the concept of aliasing and the Nyquist sampling rate. 

When dealing with dynamic real time signals, the bandwidth of the signal may be 
quite large. In order to effectively sample the signal, an arbitrary frequency limit must be 
set by the programmer. Frequency components above this limit are usually not required for 
systems analysis and may be discarded prior to sampling. One of the easiest methods to do 
this is by an analog low pass filter. Setting the low pass filter at the desired frequency limit 
allows the higher frequency components of the signal to be ignored. Because the cut oil 
frequency of a filter is defined as its half power or -3 db point, the sampling rate is usually 
somewhat greater than twice the cutoff frequency of the filter. This allows a reduced 
chance of aliasing the signal and optimum resolution in the frequency analysis. 

Another method of avoiding aliasing is to sample the signal at three to four times its 
half power frequency component. This ensures that any aliasing that takes place is in the 
region of low power signals and that it does not affect the important frequency components 
of the signal. The major drawback to this method is that by sampling at a higher rate the 
portion time signal rate is decreased. This influences resolution of the spectrum. 



x(t) 




Sampling Pulse Train tinne 




Figure 1: Sampled Analog Signal in the Time Domain 
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Figure 2: Correctly Sampled Signal in the Frequency Domain 
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Figure 3: Signal Sampled at less than the Nyquist Sampling Frequency 



When the signal is transformed into its frequency components via a discrete FFT, the 

minimum frequency resolution is inversely proportional to the time record length of the 
sampled signal. Thus if the sampling rate is F s , and the time record consist of n samples, 

then the frequency resolution will be: 



F = _J = £s. 

min n A t n 



[Eqn. 4] 



where n is the number of samples in the data set. The n A t term is the time record length 
of the rime signal and is noted by T. 
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Thus, in order to achieve the best frequency resolution, the longest possible time 
record length is desired. To do this with a fixed number of samples, the slowest sampling 
frequency possible must be used. A filtered signal in conjunction with a slow sampling 
rate will yield non-aliased, high resolution frequency information. 

Fast Fourier transform techniques assume that the time signal is periodic in nature. In 
order to meet this periodicity requirement, the time signal must be weighted, or windowed. 
This involves tapering the information by various methods without altering the frequency 
components of the signal. One of the primary windows used is the Hanning or cosine 
window. This multiplies the time signal by a sine wave of amplitude V 2/3 and a period of 
2T. The time signal is then zero in the region of t = 0 and t = T. This introduces a 
frequency component at the frequency of half of the minimum frequency resolved by the 
Fourier transform method which will not affect the results obtained by the FFT. Thus the 
periodicity requirements are met and the frequency components of the time signal are not 
altered. 

Important additional measurements performed on the mechanical signatures include 
the power spectra density function (or autospectra), the cross spectra density function, the 
system transfer function, and coherence. The input power spectra, G xx (co), is a measure of 
of the energy content of the input signal at a frequency co. Likewise, the output power 
spectra, Gyy(co), is a measure of the energy content of the output signal. Both the input 
power spectra and output power spectra are real valued functions. The cross power spectra 
G xy (co), is a complex value, 



G X y(CO) = G X yR(CO) + j G X yl(C0), 



[Eqn. 5] 



where G xy R(co) and G xy i(co) are the real and imaginary portions of the cross power 
spectra. The cross power spectra is a measure of the amount of energy transferred from the 
input point, x, to the output point, y, at a frequency co. 

The system transfer function, G(co), is the ratio of the Fourier transform of the input 



and output time signals, 

T{y(t)} Y(o» 



[Eqn. 6) 



It too is a complex number. It consist of 
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G(co) = Gr (,co) + j G kco). 



[Eqn. 7] 



where Gr(co) and Gi(co) are the real and imaginary portions of the transfer function. One 
method of calculating the transfer function uses the auto spectra and cross spectra 
information. 



|G(co)|2 = ^n 



G vv (co) 



and 



Gxx(co) 



G R (0»=^f 

Gxx(^) 

Gxx(CO) 



[Eqn. 8] 



The coherence, y X y, is a measure of the noise contamination of an output signal y(t) 
from sources other than the input signal x(t). The coherence function can be interpreted as 
the fractional portion of the output spectrum of y(t) which is linearly due to x(t) at 
frequency co [Ref. 3]. A coherence of 1 indicates perfect transmission from the input point, 
x, to the output point, y. This indicates that there is no noise contamination. A coherence 
of 0 (zero) indicates severe contamination by noise at the output point. This indicates that 
none of the signal at y is due to excitation at point x. The coherence is calculated by 



Y 2 xy(co) = 



lG X y(CO)| 2 



Gxx(^) Gyy(CO) 



[Eqn. 9] 



B. THE HILBERT TRANSFORM AND THE FREQUENCY RESPONSE 

The Hilbert transform has been identified as a means of determining the degree of 
linearity of a system's frequency response. It has been shown that accurate identification 
of a range of commonly occurring structural non-linearities such as friction, stiffness and 
power-law damping is possible in systems where the modes are well separated [Ref. 4]. 
By taking the Hilbert transform of a systems frequency response and comparing it to the 
original frequency response, the linearity of the system can be readily determined. For 
linear systems the Hilbert transform of the frequency response H(co) will be identical to the 
frequency response of the system G(co). 

The Hilbert transform of a real-valued function g(t) is defined by: 
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oo 




-oo 



[Eqn. 10] 



In the frequency domain it is defined by: 

OO 

H(co) = PV [Eqn. 11] 

J C0-^ 

-oo 

where H(co) is the Hilbert transform of the complex valued G(^) and PV is defined as the 
Cauchy Principle Value for the integral. 

G© = G R (5)+jGi«) [Eqn. 12] 

where Gr© and Gi© are the real and imaginary parts of the G©. It is important to note 
that the Hilbert transform is in the same domain as the original function. Thus the Hilbert 
transform of a real valued time function is in the time domain, and the Hilbert transform of 
a function in the frequency domain is in the frequency domain. The Hilbert transform 
H(co) is complex-valued in the frequency domain, 

H(co) = Hr(co) + j Hi(co), [Eqn. 13] 



where Hr(co) and Hi(co) are the real and imaginary parts of H(co). 

In order to avoid the singularity point at % = co, the Hilbert transform the integration 
must be divided into two portions with limits as follows: 

oo o y~ £ oo 



Hto) = pv PSr d ^ = 5 /-rr^ + /^r d? • [Eqn - 141 



-oo 



-oo 



co+ e 



Because Gr(^) is an even function and Gi(£) an odd function, 



G(-0 = G*© f 



[Eqn. 15] 
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where * denotes the complex conjugate of the complex valued function G(£), Equation 14 
must be further subdivided into its real and imaginary parts: 



H r (co) = --PV 

71 



1 Gi(ji) 
^ 2 -co 2 






oo 

Hj(co) = - — PV f ° R( ^, [Eqn. 1 6] 

it J -co- 

One method utilized to calculate the Hilbert transform involves a numerical 
integration. Due to the finite frequency resolution of a discrete Fourier transform, a 
numerical integration of the frequency response will always yield inaccuracies at the 
singularity point oo = The equations for evaluating the Hilbert transform numerically are 
given by [Ref. 5]: 



H R (coj) - K 






n 

X 

k=1 

n 

X 



G,(co k )co k Aco 

o + E. 



COk - CO j 



G R (co k ) Aco 

2 2 
COk - COj 



+ E 



[Eqn. 17] 



k-i 



where H(oo) denotes the Hilbert transform and G(co) is the system's transfer function. The 
R I 

correction terms, E j and E • , are dependent upon the system’s frequency response and the 



range of the frequency response under examination. 

A correction term is needed to correct for resonant peaks outside the range under 
examination, and others to correct for the truncation of the frequency response information 
at the limits of the frequency band under examination. Corrections for multimode systems 
are usually limited to a shift of the imaginary portion of the Hilbert transform. The 
estimate of the shift necessary to correct this is given by [Ref. 6]: 
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[Eqn. 18] 




where C Is is the correction to the imaginary portion of the Hilbert transform due to the s^ 



resonant frequency, co s . X s is the the imaginary portion of the modal constant for the s 1 * 1 
resonant frequency. An estimate of the value of X s may be evaluated from the the fact that: 



5 S , the structural damping ratio, may be estimated by using the half power bandwidth of the 
resonant frequencies outside the band under investigation. The value of Hj(co s ) may be 

taken directly from the imaginary portion of the transfer function. Thus as the resonant 
frequencies outside the band under investigation, co s , get larger, their influence of it 
correction term C Is diminishes. For structures with many resonant modes, only those 
modes near the frequency band under investigation need be used as correction terms. 

A second method of evaluating the Hilbert transform applied in this study is the use 
of the Fourier transform to perform the Hilbert transform. The integral 



is the convolution integral in the frequency domain. Since convolution in the frequency 

domain is equivalent to multiplication in the time domain, this is equivalent to taking the 
Fourier transform of G(co) and multiplying it by the Fourier transform of — — . 

7UC0 

The Fourier transform of — — is: 




[Eqn. 19] 



thus 



X s = H,(C0 S ) 8 S co. 



[Eqn. 20] 



CO 




[Eqn. 21] 



-CO 



7U CO 




t>0 



[Eqn. 22] 



The Fourier transform of G(£) is: 



r[G£)]=g(t), 



[Eqn. 23] 
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where both G(%) and g(t) are complex valued. Multiplying the Fourier transform of the 

transfer function, G(^), by the Fourier transform of— — yields: 

K CO ' 






7 "S 



if 









= 



-j g(t) t>0 
0 t = 0 . 

j g(t) t<0 



[Eqn. 24] 



Multiplying by j yields the j times the Fourier transform of the Hilbert transform: 






r { j h (co) } = j r S — p— ~ di; 

n J co -q 



= 



g(t) t>0 
0 t = 0 

-g(t) t<0 



[Eqn. 25] 



Applying the property of linearity of the Fourier transform and adding the Fourier 
transform of the transfer function G(^) to j times the Fourier transform of the Hilbert 
transform yields: 



r 2g(t) t>0 

P{G(co) + j H(co)} = g(t) t =0 , 

l 0 t<0 



and 



G(co) + j H(co) = 7 




2 

1 

0 




[Eqn. 26] 



[Eqn. 27] 



There exist many fast Fourier transform routines for complex vectors in computer 
libraries. Thus an easier method to calculate the Hilbert transform which avoids the 
numerical integration required of the Cauchy Method is to take the Fourier transform of a 
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I 2 t>0 

real signal G(co) and multiply it by{ 1 t=0 . Then by taking the inverse Fourier 

l 0 t <0 

transform of this modified step response, the transfer function G(co) and its Hilbert 
transform H(co) can be separated out [Ref 7]. Figure 4 summarizes the method used to 
calculate the Hilbert transform using the Fourier transform, Equations 21 through 27. 




Figure 4: Fast Hilbert Transform Flow Diagram 



When using the FFT to calculate the Hilbert transform a FFT routine capable of 

performing inverse Fourier transforms must be used. Because the FFT routine assumes 

convolution in the frequency domain has taken place, the Fourier transform is symmetric 
T 

about the point y, where T denotes the time record length of the input signal. The inverse 

Fourier transform also assumes that convolution has taken place. An equivalent method of 

T 

setting all negative frequency component to zero is by setting all components above y equal 

T 

to zero. All component from f>0 to y are multiplied by 2. The inverse Fourier transform 
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is then performed yielding the original function and its Hilbert transform. This method is 
faster and just as accurate as the numerical integration method for systems when applied 
properly. Correction terms are still required when using the FFT method of calculating the 
Hilbert transform. 
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III. EXPERIMENTAL PROCEDURE 



A. EQUIPMENT UTILIZED 
1. The MASSCOMP 5400 

The MASSCOMP 5400 computer used to acquire the data and process it for the 
experimental work was configured as follows: 

171 Mbyte hard disc drive 
1 300 kbyte "floppy" disc drive 
6 Mbyte Random Access Memory 

1 16 channel, 1 MHz maximum aggregate sampling rate 12 bit 
resolution Analog to Digital Converter [AD12FA] 

1 16 channel, 333 kHz maximum aggregate sampling rate 12 bit 
resolution Analog to Digital Converter [EF12M] 

5 Timing Clocks, maximum frequency 6 MHz 
1 (2) channel Digital to Analog Converter [EF12M] 

1 16 channel sample and hold temporary storage card [SHI 6 A] 

Various display and output devices 



The capabilities of the various data acquisition devices installed in the computer are 
summarized in Appendix A. 

The MASSCOMP 5400 is a Motorola 68020 based 32 bit word length system. 
It is run on the UNIX operating system. Language capabilities include "C" and 
FORTRAN. Software use to operate the system includes a graphics presentation package, 
Dynamic Signal Analysis Subroutines, Window Presentation Packages and the software 
necessary to operate the Data Acquisition devices. 

The data analysis routines, except where noted, were copied from those 
provided by the MASSCOMP Dynamic Signal Analysis Package (SP-60). These were 
originally copied from the collection of Programs for Digital Signal Processing from the 
IEEE Acoustics, Speech and Signal Processing Society. 

The programming necessary to collect and process the data was written in 
FORTRAN. All programs referred to in the remainder of this paper are documented in the 
appendices. It should be noted that the software packages for the data acquisition and 
presentation portions include subroutine calls that are unique to the MASSCOMP system 
and are written specifically to support MASSCOMP hardware. Further investigation into 
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the software documentation provided by MASSCOMP is warranted for a complete 
understanding of the data collection and presentation routines. 

In addition to the MASSCOMP specific equipment, various recording devices, 
dynamic signal analyzers and signal generators were used to provide the necessary data and 
a means of verification of the accuracy and applicability of the MASSCOMP 5400 system. 

2. Analog to Digital Converters 

Analog to Digital converters work in several different manners. While the 

method used to digitize a signal is not important, the ultimate goal of a converter is to 

compare a sampled voltage and to digitize it to the nearest voltage value within its dynamic 

range. A twelve bit A/D converter is able to separate its dynamic range into 2 12 individual 

20 

levels. Thus a converter with a 20 volt dynamic range is able to obtain a resolution of 

or 0.00488 volts. A 12 bit converter with only a 10 volt range will have twice this 
resolution. The speed of the conversion is dependent upon the method of digitization; the 
fastest devices having speed on the order of 20 nanoseconds and the slower ones of 1.2 
microseconds [Ref. 8 p. 97]. 

The transfer rate of the data to a usable state is dependent upon the 
programming of the converter. There are three methods available to the programmer for 
data acquisition. The first is reading the information directly into an array for later use in 
the program. This is the easiest and usually the quickest. If the data is read into dynamic 
memory then the converter can be run at its maximum speed or at the data bus’ speed 
whichever is less. The second method is a direct disc transfer. This is usually slower than 
the array transfer method. If the information is to be read into a hard disc for later retrieval 
and reduction then the programmer is usually limited to the transfer rate capability of the 
particular hard disc system available. One drawback of the direct to disc transfer is that the 
programmer has no access to the information until after the transfer is complete. The third 
and more complicated method is to read the data into memory buffers. The advantage of 
this method is that the data can be manipulated as the buffers are filled rather than having to 
wait until the final array or disc transfer is complete. Most micro computers allow a 
programmer to select any one of the three possibilities depending upon the application of 
the data. 

Once the data was transferred to the hard disc in either the raw or reduced form 
it could be transferred to magnetic tape or floppy disc. This allows the volume of the 
information on the hard disc to be substantially reduced after the acquisition and analysis is 
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completed. If the data is required at a later date it is available for transfer from the 
permanent recordings (floppy disc or tape). 

3. Digital to Analog Converters 

Digital to Analog converters work in somewhat the same manner of the A/D 
converter, only their task is simpler. Because the goal is to send out an analog voltage, 
there is no need to compare the desired value to the output signal. Because the D/A 
converter is essentially a resistor capacitor circuit, the settling time requirement of the RC 
circuit, and thus the D/A converter determines its speed capability. Settling times are 
dependent upon the exact circuitry and the resolution of the converter. Eight bit converters 
may have setding times of 0. 1 microseconds, while sixteen bit converters may have settling 
times of 75 microseconds [Ref. 8 p. 95]. The dynamic range of D/A converters are 
comparable to that of A/D converters. The programmer may also use the same three 
methods of transfer (array, disc and buffer) for digital to analog conversions. 

Analog to Digital converters and Digital to Analog converters require timing 
clock signals in order to determine the frequency of sampling or output. Programmable 
clock chips are used to provide voltage pulse trains for the devices. Thus the programer is 
able to set the sampling rate of an A/D device and the output speed of a D/A device. With 
clock chips capable of providing 6 MHz pulse trains, the programmer is limited only by the 
time characteristics of his conversion devices and the computers systems data bus capacity. 

4. Verification Procedure 

The process of determining whether the MASSCOMP 5400 system was 
suitable for the near real time processing of shock and vibration signals was completed in a 
series of steps. First, using analog data pre-recorded on magnetic tape from underwater 
explosions, the analog converters were tested for accuracy and speed. The MASSCOMP 
sampled, recorded, analyzed and displayed the information. Comparisons to the data 
analyzed by other means were used to verify the accuracy of the acquired data. Once 
validation of sampling speed and analysis of selected data sets was accomplished all data 
sets were transferred to the hard disc for archival purposes. 

The second method of systems verification was completed by impact testing of 
an aluminum plate. Data analysis was expanded to include power spectral analysis and 
coherence measurements. Again the results were compared to those obtained by an 
separate HP 3562A Dynamic Signal Analyzer in order to verify the accuracy of the 
MASSCOMP system. 
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5. A pplication of the MASSCOMP 5400 System 

Once the system was determined to operate satisfactorily, it was used on the 
vibrational analysis of a linear response system. In addition to the frequency response, 
power spectral density and coherence measurements, the Hilbert transform was used to 
examine the degree of linearity of the frequency response. This enabled a rapid 
identification of any nonlinearities in the system transfer function. The two methods of 
computing the Hilbert transform, numerical integration and use of the Fourier transform, 
were compared for accuracy. Once the Hilbert transform method was chosen, the 
frequency response characteristics of a non-linear system were evaluated to ensure the 
Hilbert transform identified the non-linearities. 

B. SHOCK AND VIBRATION ANALYSIS 

1. Analysis of Prerecorded Underwater Explosion Data 

The analysis of prerecorded underwater explosion data involved the transfer of 
pressure and strain gage records from a Honeywell Model 101 sixteen channel FM tape 
recorder to the MASSCOMP system. Data was recorded at a tape speed of 120 inches per 
second. The AD12FA Analog To Digital converter card was used for the transfer of data. 
Figure 5 shows the schematic diagram of the connections used for the data transfer. 

The bandwidth capability of the tape recorder was 250 kHz. The bandwidth of 
the pressure transducer data was 250 kHz, while the strain gage data was only 10 kHz. 
Location of data on the tape was noted by a trigger signal on one of the sixteen data 
channels. During the explosive testing procedure, an electrical wire with a DC voltage 
applied to it was wrapped around the explosive charge. As the charge was detonated, the 
cable was severed and the voltage sensed at the recorder dropped to zero. Monitoring the 
trigger channel for this drop in the voltage enabled the program to start the clock module. 
When the clock began sending pulses to the Analog To Digital converter, information was 
transferred off the tape into the data storage arrays in the computer program. 

The Analog to Digital Converter section of the EF12M multifunction card is 
limited to a bandwidth of only 160 kHz. This required the use of the AD12FA Analog To 
Digital converter. This allowed sampling of the recorded data at a maximum frequency of 
1.0 MHz, well beyond the bandwidth of the original signal. Figure 6 is the pressure 
transducer data set used as a reference to compare the data recorded by the A/D system. 
This was recorded by a HP5451C computer system at a sampling rate of 500 kHz. Data 
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was acquired on the MASSCOMP system at various sampling rates as shown in Table 1 
and then compared to the reference record. The program utilized to acquire the data is 
shown in Appendix B. 

TABLE 1: 

SUMMARY OF SAMPLING RATES USED TO EXAMINE UNDERWATER 

EXPLOSION DATA 



Tape Speed (ips) 


Sampling Rate 


120 


1 MHz 


120 


750 kHz 


120 


500 kHz 


120 


250 kHz 




Figure 5: Schematic Diagram of Connections Used for Underwater Explosion Data 

Transfer 
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Figure 6: Reference Pressure Record Used as the Standard 

Figures 7 through 10 are the data recordings from the same pressure pulse 
recorded at the sampling rates listed above. Comparisons of these Figures to the base line 
record revealed that the 500 kHz sampling rate matched the reference ideally. This was as 
predicted by the Nyquist sampling theorem. Recordings at rates less than the Nyquist 
frequency missed portions of the data that was available on the tape. Rates faster than the 
Nyquist frequency exceeded the bandwidth of the tape recorder and the reference record 
and subsequently appeared "noisy". Because there was no method to determine if this 
noise was actually recorded information or noise produced by the tape recording and play 
back process, the sampling rate for all further analysis of the prerecorded underwater 
explosion data was fixed at 500 kHz. The strain gage data was then acquired from the 
Honeywell tape recorder and transferred to the MASSCOMP hard disc. 

The data was transferred to the the MASSCOMP system one track at a time. In 
order to utilize the system, as it is presently configured, for real time data acquisition in the 
field the sixteen channel sample and hold card (SH16F) and the single function Analog to 
Digital Converter (AD12FA) must be used. The only limitation is the data transfer bus, 
which is limited to 2 Megabytes per second This allows one million sixteen bit words per 
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second to be gathered in real time. Table 2 summarizes the bandwidth limitations based 
upon the number of sampled channels. 

TABLE 2: BANDWIDTH LIMITATIONS FOR MULTIPLE DATA CHANNELS 



Number of Channels Sampled 


Bandwidth Limit 


1 


500 kHz 


2 


250 kHz 


4 


125 kHz 


8 


62.5 kHz 


16 


31.5 kHz 



Thus it would be impossible to record more than two pressure pulses in real 
time at a 500 kHz sampling rate. It would be possible to record up to sixteen strain gage 
traces at a 20 kHz sampling rate on the AD12FA A/D. The present system configuration 
limits the data acquisition to one type of data, either pressure or strain gage. The single 
sample and hold card can only be clocked at one frequency, limiting the choice to one type 
of data. A second sample and hold card, connected to the EF12M multifunction module, 
would allow the simultaneous sampling of both pressure and strain gage information. 

In addition, to utilize the system for live data recording the program listed in 
Appendix B would have to be modified to include the use of the sample and hold card. 
This would necessitate a modification to the identification of the channels to be sampled and 
the clocking method used. Due to the increase in the volume of data, the storage arrays 
would have to be modified to allow up to sixteen channels of data. This would create the 
need for large arrays (16 [one for each channel] by 20,000 [one second of data at 20 kHz 
sampling rate]) for the raw data. Another method would be to read the information directly 
to the disc and then process it after the data is on the disc. This would still allow 
processing of the data in the field and allow analysis within minutes of the explosion. 
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Figure 7: Pressure Recording at 1 MHz Sampling Rate 
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Figure 8: Pressure Recording at 750 kHz Sampling Rate 
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Figure 9: Pressure Recording at 500 kHz Sampling Rate 
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Figure 10: Pressure Recording at 250 kHz Sampling Rate 
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2. Impact Testing of an Aluminum Plate 

Using the equipment schematic diagram shown in Figure 1 1, the digital signal 
analysis portion of the data acquisition system was tested. A modally tuned impact hammer 
with a bandwidth of 1,200 Hz was used to strike the aluminium plate. A force transducer 
located in the hammer provided a voltage signal to one channel of the A/D converter. An 
accelerometer located on the plate provided the plate acceleration data to a second channel of 
the A/D converter. 

The data was sampled at a rate of 5,000 Hz per channel providing a 2,500 Hz 
bandwidth for the sampled signal. Because the energy level of the input signal was 
negligible beyond 2,500 Hz, a low pass filter was not required. Data blocks of 1024 
samples were used yielding a frequency resolution of 4.88 Hz. This ensured that an 
adequate number of data points were in the frequency band width of the hammer. The 
program utilized to acquire the data is shown in Appendix C. 

The time signals for the hammer and plate were recorded and averaged in the 
frequency domain. The frequency information was derived utilizing the FFT routine 
FFT842 [Ref. 9], while the power spectral information and coherence measurements were 
derived by utilizing the CCSE routine [Ref. 10]. Both of these routines were modified to 
suit the particular characteristics of the impact testing procedure and the microcomputing 
system. 

The dynamic signal analysis was first performed on a single data record. The 
coherence, transfer function and frequency components of the signals were computed. As 
expected the coherence for a single impact was 1 throughout the frequency band under 
examination. The identical test was recorded using a HP3562A Dynamic Signal Analyzer 
for comparison. The analysis performed utilizing the super microcomputer yielded 
comparable results. Figures 12 through 14 display the recorded data for this single impact 
test from both sources. 
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Figure 11: Experimental Set Up for Impact Testing of Aluminum Plate 



The impact testing was then averaged to determine the influence of averaging on 
the transfer function. Figures 15 through 18 are the data recordings for sixteen impacts. 
As expected the frequency components of both the hammer and accelerometer data 
smoothed considerably when averaged. The coherence decreased in areas of noise 
contamination and non-linearities. 

3. Steady State Vibration Response Measurement 

Steady state vibrations were measured using the apparatus shown in the 
schematic diagram in Figure 19. An aluminum beam was mounted so that it could be 
excited by a shaker. The input acceleration was measured at the mounting jaws and the 
output acceleration was measured at the tip of the beam. The beam was excited using a 
random noise excitation of 1 kHz and 500 Hz bandwidths. The system transfer function 
was calculated, and the Hilbert transform of the real and imaginary portions of the 
frequency response was performed to determine if the system behaved in a linear manner. 
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Figure 12: Input Power Spectra for Single Impact 
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Figure 13: Output Power Spectra for Single Impact 
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Figure 14: Transfer Function for Single Impact 
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Figure 15: Input Power Spectra for Multiple Impact 
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Figure 16: Output Power Spectra for Multiple Impact 
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Figure 17: Transfer Function for Multiple Impact 
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Figure 18: Coherence Measurement for Multiple Impact 

A low pass filter ensured that no frequency components above the range of 
investigation contaminated the transfer function. The sampling rate was such that a 
frequency resolution of 1.3 Hz was obtained. For steady state vibration response 
measurements, a modified version of the program used for impact testing was used. The 
program is listed in Appendix D. Because the same data reduction and analysis was 
performed on the steady state vibrations testing as was done on the impact testing 
information the same reduction routines were used. 

Figure 20 is the input power spectra for a random noise excitation at a 1 kHz 
bandwidth centered at 500 Hz. Figures 21 and 22 are the output power spectra and 
coherence measurements for the same excitation. The transfer function for this excitation is 
shown in Figure 23. The use of random noise produced a low coherence measurement in 
the region above 200 Hz. In order to improve the coherence measurements and smooth the 
transfer function in the low frequency area, a random noise excitation of only 500 Hz 
bandwidth starting at zero was used. The results of using this excitation are shown in 
Figures 24 through 29. 

An attempt to utilize a swept sine wave as the excitation source was made. Two 
methods, one using a HP 3562A Dynamic Signal Analyzer and a second using a voltage 
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controlled oscillator driven by the digital to analog converter of the EF12M multifunction 
module, were used to sweep through a frequency band. Using the Hewlett Packard 
dynamic signal analyzer was not effective. As the dynamic signal analyzer swept through 
the frequency bandwidth, the MASSCOMP would sample the responses. However, due to 
the length of time required by the MASSCOMP to process the data, the MASSCOMP 
would miss some frequencies as the HP Dynamic Signal Analyzer swept through them. 
The method of using a voltage controlled oscillator was more effective; however, the 
processing time was such that only a 200 Hz bandwidth could be swept without skipping 
some frequencies. This required approximately twenty-five minutes to sweep this 
bandwidth, while with a Dynamic Signal Analyzer the entire bandwidth from 0 to 1 kHz 
could be sampled in the same amount of time. This necessitated the use of random noise as 
the only excitation source of the shaker. 

The results of the steady state analysis were compared to the HP 3562A 
Dynamic signal Analyzer for accuracy. Figures 30 through 34 are the results of those test. 
When compared to the results gained via the microcomputer, the differences are quite 
noticeable. Of particular importance is the difference in magnitude in the input power 
spectra, the output power spectra and the cross spectra. While the shape of the 
measurement results are nearly identical, the magnitudes are different. This leads to an 
incorrect approximation of the real and imaginary portions of the transfer functions. This is 
evident in Figures 35 through 38. Figures 35 and 36 are the real and imaginary portions of 
the transfer function measured with the microcomputing system. The identical 
measurements using the dynamic signal analyzer are shown in Figures 37 and 38. The 
difference in magnitude is immediately apparent. 
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Figure 19: Steady State Vibration Measurement Schematic Diagram (upper) and 

Cantilevered 

Beam used in Experimental Procedure Gower). 
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Figure 20: Input Power Spectra for 1 kHz Random Noise 
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Figure 21: Output Power Spectra for 1 kHz Random Noise Excitation 
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Figure 22: Coherence Measurement for 1 kHz Random Noise Excitation 
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Figure 23: Transfer Function for 1 kHz Random Noise Excitation 



39 



-15 



Gxx (f) 




FREQ (Hz) 

RtfXJO-500 



Figure 24: Input Power Spectra for 500 Hz Random Noise Excitation 
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Figure 25: Output Power Spectra for 500 Hz Random Noise Excitation 
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Figure 26: Coherence Measurement for 500 Hz Random Noise Excitation 
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Figure 27: Transfer Function for 500 Hz Random Noise Excitation. 
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Figure 28: Real Portion of Cross Power Spectra for 500 Hz Random Noise Excitation 
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Figure 29: Imaginary Portion of Cross Power Spectra for 500 Hz Random Noise 

Excitation. 
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Figure 30: Input Power Spectra for 500 Hz Excitation Measured with HP 3562A 




Figure 31: Output Power Spectra for 500 Hz Random Excitation Measured with HP 

3562A 
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Figure 32: Transfer Function for 500 Hz Random Noise Excitation Measured with HP 

3562A 
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Figure 33: Real Portion of Cross Power Spectra for 500 Hz Random Excitation Measured with HP 3562A 
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Figure 34: Imaginary Portion of Cross Power Spectra for 500 Hz Random Excitation Measured with 

HP 3562A 
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Figure 35: Real Portion of Transfer Function for 500 Hz Excitation 
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Figure 36: Imaginary Portion of Transfer Function for 500 Hz Excitation 
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Figure 37: Real Portion of Transfer Function for 500 Hz Excitation Measured with HP 

3562A 
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Figure 38: Imaginary Portion of Transfer Function for 500 Hz Excitation Measured with 

HP 3562A 
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IV. APPLICATION OF THE HILBERT TR ANSFORM 



The aluminum beam used for steady state vibrations analysis was also used to 
evaluate the feasibility of including the Hilbert transform as a method of identifying system 
irregularities. There are two resonant modes, one at 137 Hz and the second at 891 Hz, 
shown previously in Figure 23. Using the half power bandwidth method to estimate the 
viscous damping ratio at the second resonant frequency: 




- co 1 ^ 

^ 0)5 J 



[Eqn. 28] 



where co 2 and cOj are the half power frequencies for that resonant point and co s is the 

resonant frequency. Using the information for the second resonant frequency from Figure 
39, the viscous damping ratio was determined to be 0.00180. This damping ratio was 

then used to estimate the correction terms necessary for the Hilbert transform. Recalling 
that the structural damping ratio is proportional to the viscous damping ratio, 

5 S = f . [Eqn. 29] 

and substituting into Equation 20, 

X s = H,(o) s ) 8 S 0 ) s =H I (to s )fto s , [Eqn. 30] 



Substituting Equation 30 into Equation 18, the correction term for the imaginary portion of 
the Hilbert transform is then: 



C Is (CD)=j 



CO X s 

2 2 
CO s - CO 



= j 



CO H I (co s ) co s 
2 (COs - CO 2 ) 



CO Hj(co ) (co 2 - CO,) 

j 2 2 

4 (co s - co 2 ) 



[Eqn. 31] 



Using the damping information from the second resonant frequency, the correction terms 
were less than 1 0' 5 and were considered negligible and were not included. 
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The real and imaginary portions of the original transfer function and the Hilbert 
transform of the transfer function are shown in Figures 40 through 44. The numerical 
integration method of performing the Hilbert transform was not accurate in the area of the 
first resonant frequency. Figures 45 and 46 show that the Hilbert transform of the transfer 
function was not significantly different than the original transfer function. The numerical 
integration method will always yield incorrect answers for systems of low damping at 
resonance. 

The FFT method of performing the Hilbert transform was also inaccurate. Due to the 
transform into the time domain and subsequent truncation, there appeared to be high 
frequency components introduced to the transfer function. This is evident in Figures 47 
and 48. However, the shape of both the real and imaginary portions of the Hilbert 
transform are identical to the original transfer function. 

When using the frequency response of a model of a linear system the Hilbert 
transform was able to indicate that the system behaved in a linear fashion. However, due 
to the inaccuracy in measuring the cross power spectra the transfer functions used in 
performing the Hilbert transform are inaccurate. Thus in order to ensure that the Hilbert 
transform in yielding correct informadon, the original transfer function must be correct; 
otherwise, the Hilbert transform will indicate a system that is in fact a linear system as a 
non-linear system. 
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Figure 39: Second Resonant Frequency Identified using 1 kHz Random Noise 

Excitation 
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Figure 40: Real Portion of the Transfer Function of the First Resonant Frequency 
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Figure 41: Imaginary Portion of the Transfer Function of the First Resonant Frequency 
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Figure 42: Hilbert Transform of the Real Portion of the Transfer Function using Numerical Integration 
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Figure 43: Hilbert Transform of the Imaginary Portion of Transfer Function using Numerical Integration 
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Figure 44: Hilbert Transform of the Real Portion of the Transfer Function using FFT Method 
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Figure 45: Hilbert Transform of the Imaginary Portion of the Transfer Function using FFT 

Method 
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V. CONCLUSIONS 



A. UNDERWATER EXPLOSION TESTING 

1 . The super microcomputer is able to record up to two channels of shock pressure with 
250 kHz bandwidth or sixteen channels of strain information with 10 kHz bandwidth 
directly to data arrays. It is possible to reduce and display the data prior to storage to 
a recording media. 

2. In the present configuration of only one sample and hold device for the data 
acquisition hardware only one type of data, shock pressure or strain, can be recorded. 
The addition of a second sample and hold card would allow two different data types 
to be recorded and analyzed. 

B . IMPACT TESTING OF AN ALUMINUM PLATE 

1 . The data acquisition capabilities, tied to the Digital Signal Processing software, allow 
for a complete dynamic signal analysis of a structure using a microcomputer. 

2. Sampling rates of 32 kHz are sufficient for multi-channel low frequency vibrations 
analysis. However, the requirement to use an analog filter prior to sampling requires 
the use of a filter for each data channel monitored. This would require a rack of 
filters and reduce the ease of movement of the system. Setting up of a vibrations 
analysis station on a permanent basis would be required. 

C. STEADY STATE VIBRATION ANALYSIS 

1 . The data acquisition techniques are sound; however, the microcomputer data 
acquisition equipment yields information that is significantly different in magnitude 
when compared to a proven measurement device. This is due to the 12 bit resolution 
over the 10 volt dynamic range of the analog to digital converter. 

2. The coherence measurements using an external random noise source indicate that this 
is not the optimum method of performing response measurements. Attempts to 
improve the coherence measurements by using a voltage controlled oscillator 
controlled by the computer as a sinusoidal excitation source were unsuccessful. 

3. The identification of non-linear response using the Hilbert transform is limited to 
systems with a few widely spaced resonant modes. It requires data reduction to be 
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performed by the user prior to attempting the identification of non-linearities. 
Because of the poor voltage resolution of the analog to digital converters, the system 
transfer information used as an input for the Hilbert transform is in error. Thus the 
reliability of the Hilbert transform for use determination of the linearity of the 
frequency response is poor. 

4. The present data presentation in graphical format capability is rudimentary and 
requires improvement. The capability to include cursor input into the program for 
axes scaling and data zooming and recording would greatly improve the analysis 
capabilities available to the user. The capabilities of using a "mouse" or a digital 
graphics pad as an control device to delineate points of interest for further reduction 
would greatly enhance the capability to process information. 
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VI. RECOMMENDATIONS 



A. UNDERWATER EXPLOSION TESTING 

1 . Include the use of a second sample and hold card. This would allow simultaneous 
recording of a single high frequency shock pressure signal and sixteen low frequency 
strain gage signals. 

B. IMPACT TESTING 

1 . Improve the presentation software so that cursor input via a "mouse" can be used to 
select areas of interest for further analysis. 

C. STEADY STATE VIBRATION ANALYSIS 

1 . Explore the possibility of using a proved data acquisition device, the HP 3562A, in 
conjunction with the programming capability of the MASSCOMP system. The 
EF12M Multifunction Module has an IEEE 488 interface capability. The HP 3562A 
may be able to acquire the data and transfer it to the MASSCOMP for exploration of 
the degree of linearity of the transfer function by use of the Hilbert transform. 

2. Reduce the processing time requirements for large data arrays. FORTRAN 
programming is not the best method to process information in real time applications. 
Transferring the code to assembly language or machine code would greatly reduce the 
time requirements of data analysis and allow faster structure analysis. 

3. Improve the presentation software as noted above for impact testing. 

4. Develop a method of swept sine wave excitation of the structure using the 
microcomputer as the controller. The digital to analog converter allows for two 
channels of control information to be sent from the computer, yet it must operate 
independent of the analog To digital converter. 

5. Add a signal conditioner to boost the level of the input signals to the analog to digital 
converter. This would result in more accurate calculation of the system transfer 
function and increase the reliability of the results of the Hilbert transform 
measurements. 
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APPENDIX A 



MASS COMP 5400 HARDWARE CAPABILITIES 



The capabilities of the data acquisition devices installed on the MASSCOMP 5400 and 
the microcomputer itself are summarized below: 

MASSCOMP 5400 Super Microcomputing system 
Motorola 68020 CPU module 
6 Mbyte of system memory 
Data Acquisition Hardware 

AD12FA Analog to Digital Converter 
EF12M Multifunction Module 
SH16F Sample and Hold Card 
Two Display Terminals 
HP plotter 
Dot matrix printer 



AD12FA Analog to Digital Converter 
[Path identified as /dev/dacpO/adfl] 

16 channel input single ended, 8 differential [channels 0-15] 

Programmable-gain amplifier (with gains selectable by software of 1,2,4, or 8) 
Sample and Hold circuit for connection to the SH-16F 
Maximum aggregate sampling rate: 1 Mhz 

Resolution: 12 bits (zero padding for total sample length of 16 bits) 

Input Voltages: 

0 to +10 V 
-5 to + 5 V 

EF12M Multifunction Module 

Analog to Digital Converter Section 
[Path identified as /dev/dacpO/adfO] 

16 channel input single ended, 8 differential 

Programmable- gain amplifier (with gains selectable by software of 1,2,4, or 8) 

Maximum aggregate sampling rate: 333 kHz 

Resolution: 12 bits (zero padding for total sample length of 16 bits) 

Input Voltages: 

-lOto+lOV 
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Digital to Analog Converter Section 
[ Path identified as /dev/dacpO/dafO] 

2 channels 

Maximum clock rate per channel: 333 kHz 

Resolution: 12 bits (zero padding for total sample length of 16 bits) 

Output Voltages: -10 to +10 V 
Parallel I/O Section 

16 Channels Input and Output utilizing hardware handshake for transmission 
of data 

Clock Section 

[ Path identified as /dev/dacpO/efclkn (n identifies clock 0-5)] 

5 Counters with output, source input and gate input connections 
Frequency capability ~0 Hz to 3 MHz 

SH-16F Sample and Hold Module 

[Path /dev/dacpO/adfl, channels 16 through 31 on the AD12FA card] 

16 Single-ended Input Channels or, 

8 Differential Channels 

Maximum of 3 SH-16F modules may be connected to a AD12FA or EF12M to 
give the capability of sampling a maximum of 48 channels with a single A/D 
converter. The maximum aggregate sampling rate must never exceed the 
capability of the "host" A/D converter. 
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APPENDIX B 



UNDERWATER EXPLOSION DATA TRANSFER PROGRAM 

This appendix contains the program used to collect the data of the underwater explosion 
recordings from the FM tapes. All routines were written by the author. For further 
explanation of the data acquisition specific routines the reader is encouraged to refer to the 
MASSCOMP documentation " Data Acquisition User's Manual." Figure 46 is the program 
flow diagram. Sections of the program follow the general outline presented in the flow 
diagram. Specific subroutines for plotting the information are not shown in this appendix. 



* THIS PROGRAM IS THE MAIN PROGRAM FOR DATA ACQUISITION 

* OF DATA FROM STORED TAPES OF UNDERWATER EXPLOSIONS 

* VERSION 1 

* LT R. A. SHAFER. USN 

* NAVAL POSTGRADUATE SCHOOL 

* 28 APRIL 1987 



* * * * * * * * * * * * * * * * * * * * * * * * * * * * 3ft * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

* SET UP VARIABLES 

* MYBUFSZ IS THE BUFFER SIZE FOR TEMPORARY DATA STORAGE OF A/D DATA 

* FACTOR IS A CONVERSION FACTOR FOR BUFFER MANAGMENT 



N BUFFS 
NCLKS 
NEARFREQ 
RDONLY 
STARTLO 
TIMEOUT 
FREQ 
WIDTH 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

REAL 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 



THE NUMBER OF BUFFERS FOR DATA COLLECTION 
THE NUMBER OF CLOCKS USED FOR A/D CONTROL 
PARAMETER FOR CLOCK CONTROL 
SET A/D CONVERSION TO READ ONLY 

SET CLOCK WAVEFORM IN LOW (ZERO VOLTAGE) POSITION 
TIME LIMIT FOR DATA ACQUISITION 
FREQUENCY OF SAMPLING CLOCK 

WIDTH OF SAMPLE PULSE (0.5 MICROSEC MIN) 

MYBUFSZ 
FACTOR 
NBUFFS 

NCLKS. NEARFREQ, RDONLY 
RDWRT, SLICE, SQUARE 
STARTLO, TIMEOUT 
FREQ, WIDTH 
( MYBUFSZ = 16384 ) 

( FACTOR = 2) 

( NBUFFS = 3 ) 

( NCLKS = 1 ) 

( NEARFREQ = 0 ) 

( RDONLY = 1 ) 

( RDWRT = 0 ) 

( SLICE = 10 ) 



59 



PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 



( SQUARE = 4 ) 

( STARTLO = 0 ) 
( TIME0000 ) 

( WIDTH = 0.0 ) 



SET UP INPUT BUFFER LIMITS 

ENSURE BUFFER IN EVEN WORD BOUNDRY IN STACK 

INTEGERS DATABLK(98304) 

INTEGER DUMMY 
EQUIVALENCE (DATABLK, DUMMY) 

DECLARE SERVICE ROUTINE AS EXTERNAL PROGRAM 
EXTERNAL AD SERVICE 

SET UP VARIABLES FOR A/D CONTROL AND SETTINGS 

ADPATH MARKER FOR PATH TO A/D CONVERTER 
CLKPATH MARKER FOR PATH TO CONTROL CLOCK 
NCHAN NUMBER OFA CHANNELS TO BE SAMPLED 
FCHAN FIRST CHANNEL TO BE SAMPLED 
SCHAN SPACING OF CHANNELS IF MORE THAN ONE SAMPLED 
GAIN GAIN SETTING OF A/D CONVERTER 

INDEX ARRAY POINTER FOR BUFFER MANAGEMENT 
POWER NUMBER TO RAISE 2 TO THE POWER OF FOR SAMPLE BLOCK SIZE 
TYME ARRAY FOR TIME VALUES 

RECLEN TIME RECORD LENGTH OF DATA SET 
RFREQ ACTUAL FREQUENCY CLOCK SET TO 

R WIDTH ACTUAL WIDTH OF SAMPLE PULSE 

SAMPLES ARRAY FOR TIME RECORD DATA 
LOW LOWER LIMIT FOR DATA STORAGE ON HARD DISK 

HIGH UPPER LIMIT FOR DATA STORAGE ON HARD DISK 



INTEGER 


ADPATH 


INTEGER 


CLKPATH 


INTEGER 


STATWDS(2) 


INTEGER 


I 


INTEGER 


NCHAN 


INTEGER 


FCHAN 


INTEGER 


SCHAN 


INTEGER 


GAIN 


INTEGER 


INDEX 


INTEGER 


POWER 


REAL 


TYME(16384) 


REAL 


RECLEN 


REAL 


RFREQ 


REAL 


RWIDTH 


REAL 


SAMPLES(16384) 


REAL 


LOW 


REAL 


HIGH 


REAL 


TIMSCL 


REAL 


SUM 



CHARACTER ANSWER 



COMMON SAMPLES, ADPATH, CLKPATH, INDEX, DATABLK, NITEMS 



* WRITE HEADER ON SCREEN 

* 

WRITE (6,9999) 

9999 FORMAT (47H1 UNDERWATER EXPLOSIONS DATA ACQUISITION PROGRAM) 

* 

* GET INPUT PARAMETERS FROM SCREEN 

* 

PRINT *, "INPUT NUMBER OF CHANNELS TO BE SAMPLED" 

READ *, NCHAN 

PRINT INPUT FIRST CHANNEL TO BE SAMPLED" 

READ * FCHAN 

PRINT INPUT SPACING OF CHANNELS ( 0 FOR NONE ) " 

READ *, SCHAN 

PRINT *, "INPUT GAIN FACTOR ( 0,1, 2,3 = *1, *2, * 4 , *8 )" 

READ *, GAIN 

* 

* ENSURE ALL DACP DEVICES ARE CLOSED 

* 

CALL MRCLOSALL 

* 

* OPEN A/D CONVERTER PATH 

* SET FIRST CHANNEL, SPACING AND GAIN ON A/D CONVERTER 



ADPATH = -1 

CALL MROPEN( ADPATH, "/dev/dacpO/adfO", RDONLY) 

CALL MRADMOD(ADPATH, 0, 0) 

CALL MRADINC(ADPATH, FCHAN .NCHAN, SCHAN, GAIN) 

* 

* IDENTIFY BUFFERS FOR DATA COLLECTION AND RELEASE THEM 

* TO THE EMPTY BUFFER QUEUE OF A/D CONVERTER 

* 

DO 200 I = 1, 81921, MYBUFSZ 

CALL MRBUFID (ADPATH, DATABLK(I), MYBUFSZ*2) 

CALL MRBUFREL (ADPATH, DATABLKQ) 

200 CONTINUE 
* 

* OPEN CLOCK PATH FOR A/D CONVERTER CONTROL CLOCK 

* AND INTERNALLY CONNECT TO THE A/D CONVERTER 

* 

CLKPATH = -1 

CALL MROPEN (CLKPATH, 7dev/dacp0/efclk2", RDWRT) 

CALL MR WIRE (ADPATH, 0) 

* 

* GET CLOCK SETTINGS FROM SCREEN 

* 

PRINT *, "INPUT DESIRED FREQUENCY IN HZ" 

READ *, FREQ 

* 

* DISARM CLOCK TO ENSURE IT IS NOT RUNNING 

* AND SET UP CLOCK PARAMETERS 

* 

CALL MRCLKDIS ( 1, CLKPATH ) 
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CALL MRCLK1GATED ( CLKPATH, NEARFRQ, FREQ, RFREQ, 1, 1.5, R WIDTH, 
STARTLO, 5) 

* 

* GET NUMBER OF DATA POINTS TO BE RECORDED 

* AND TIME SCALING FACTOR 

PRINT *, "INPUT POWER OF TWO FOR NUMBER OF SAMPLES ( MAX = 14 )" 
READ *, POWER 
NITEMS = 2** POWER 

PRINT *, "NUMBER OF SAMPLES IS "NITEMS 

RECLEN = NITEMS / RFREQ 

PRINT *, "ENTER TIME SCALING FACTOR" 

READ *,TIMSCL 

PRINT *, "SIMULATED SAMPLING FREQUENCY IS -> ",RFREQ*TIMSCL 
PRINT ‘."TIME RECORD LENGTH IS ".RECLEN/TIMSCL 
PRINT *, "START TAPE NOW- 

PRINT ‘."THEN ENTER "y" AND HIT RETURN WHEN TAPE IS READY" 

READ *, ANSWER 

* 

* ARM CLOCK AND WAIT FOR TRIGGER TO DROP TO ZERO 

* 

PRINT * ."ARMING CLOCK" 

CALL MRCLKARM ( 1, CLKPATH ) 

PRINT ‘."WAITING FOR TRIGGER PULSE" 

* 

* COLLECT DATA SET AND SEND BUFFER TO A/D SERVICE ROUTINE 

* 

CALL MRXINQ ( ADPATH, NITEMS, NITEMS, ADSERVICE ) 

CALL MREVWT ( ADPATH, STATWDS, 30000 ) 

* 

* COMPUTE AVERAGE OF FIRST 1200 SAMPLES TO REMOVE ANY ZERO OFFSET 

* 

PRINT *, "STOPPED READING" 

DO 250 1= 1, 1200 
SUM= SUM + SAMPLES(I) 

250 CONTINUE 

SUM = SUM/ 1200 

* 

* FILL TIME ARRAY AND REMOVE ZERO OFFSET FROM RECORD 

* 

DO 300 I = 1, NITEMS 
TYME(I)=(I- 1)/(RFREQ*TIMSCL) 

SAMPLES (I) = SAMPLES (I) - SUM 
300 CONTINUE 

9990 FORMAT (IX, G20.8, 10X , G20.8) 

* 

* CALL PLOTTING ROUTINE FOR DISPLAY OF DATA 

* 

CALL PLOTl( NITEMS, SAMPLES, TYME, " Time History " , "Voltage" , "Sccon 

ds'V'iimhis.g" ) 

* 

* WRITE DATA TO HARD DISK 

* 

* PRINT *, " INPUT LIMITS OF HARDDISK DATA RECORD " 

* READ *, LOW, HIGH 

* DO 350 I = 1, NITEMS 
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* EF ( TYME(I) .GE. LOW .ME(I) .LE. HIGH) THEN 

* WRITE (7,9990) TYME(I), SAMPLES® 

* END IF 
*350 CONTINUE 

STOP 

END 

SUBROUTINE ADSERVICE 
INTEGER K 
INTEGER INDEX 
INTEGER POWER 
INTEGER ADPATH 

INTEGER NITEMS 
INTEGER* 2 DATABLK(98304) 

INTEGER GETINDEX 
PARAMETER (GETINDEX = 3) 

REAL DUMMY 

REAL SAMPLES( 16384) 

COMMON SAMPLES, ADPATH, CLKPATH, INDEX, DATABLK, NITEMS 

IF (.EQ. 0) THEN 

NITEMS = 16384 
END IF 

* 

* STOP CLOCK, GET BUFFER FROM DONE BUFFER QUEUE 

* CORRECT DATA BY SCALING FACTOR (20/2*12) 

* FILL SAMPLE ARRAY WITH TIME RECORD AND RELEASE BUFFER TO 

* READY BUFFER QUEUE 

* 

CALL MRCLKDIS( 1, CLKPATH ) 

CALL MRBU GETINDEX, INDEX) 

TEMP = 1 / 204.8 
DO 300 K = 1, NITEMS 
DUMMY = FLOAT® AT ABLK(K)) * TEMP 
SAMPLES(K) = DUMMY 
300 CONTINUE 

CALL MRBUFREL (ADPATH, DATABLK(INDEX)) 

RETURN 

END 
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Figure 46: Underwater Explosion Data Transfer Program 
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APPENDIX C 



IMPACT TESTING DATA ACQUISITION AND REDUCTION PROGRAM 

This section contains the program listing used to acquire the data from the aluminum 
plate. All routines were written by the author except where noted. For further explanation 
of the data acquisition specific routines the reader is encouraged to refer to the 
MASSCOMP documentation "Data Acquisition User's Manual". Figure 47 is the program 
flow diagram. Sections of the program follow the general outline presented in the flow 
diagram. Specific subroutines for plotting the information are not shown in this appendix. 



* THIS PROGRAM IS THE MAIN PROGRAM FOR DATA ACQUISITION 

* AND REDUCTION OF IMPACT TESTING OF AN ALUMINUM PLATE 

* 

* R. A. SHAFER, LT, USN 

* NAVAL POSTGRADUATE SCHOOL 

* MONTEREY, CA 93940 

* 27 MARCH 1987 



* SET UP CONSTANTS AND VARIABLES 

* 



INTEGER 
INTEGER 
INTEGER 
INTEGER 
REAL 
REAL 
INTEGER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
INTEGERS 
INTEGER 
EQUIVALENCE 
EXTERNAL 
INTEGER 



MYBUFSZ, NBUFFS 
NEARFREQ, RDONLY 
RDWRT, SQUARE 
STARTLO, TIMEOUT 
FREQ, WIDTH 
PI 

OUTFILE 

( NBUFFS = 10 ) 

( NEARFREQ = 0 ) 

( RDONLY = 1 ) 

( RDWRT = 0 ) 

( SQUARE = 4 ) 

( STARTLO = 0 ) 

( TIMEOUT = 30000 ) 

( WIDTH = 0.0 ) 
DATABLK(100000) 
DUMMY 

(DATABLK.DUMMY) 
ADSERVICE 
ADPATH 
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INTEGER CLKPATH(2) 

INTEGER STATWDS(2) 

INTEGER I 

INTEGER NCHAN 

INTEGER FCHAN 

INTEGER SCHAN 

INTEGER GAIN 

INTEGER INDEX 

INTEGER 
REAL TPI 

REAL WEGHT(513) 

REAL TYME(1024) 

REAL RECLEN 

REAL RFREQ 

REAL R WIDTH 

REAL XTR(1024), Y1TR(1024) 

REAL XTI(1024) 

REAL Y1TI(1024) 

REAL XX(1024) 

REAL YY 1(1024) 

REAL XFR ( 1 024 ) ,XFT( 1 024) 

REAL Y1FR(1024),Y1FI(1024) 

REAL GXX(513) 

REAL GYY(513) 

REAL GXYRE(513) 

REAL GXYIM(513) 

INTEGER NOS AMP 

REAL DELTAFRQ(1024) 

REAL DFREQ 

REAL HFR(513),HFT(513) 

REAL HF(513),HMR(513) 

REAL HMI(513),HM(513) 

REAL AVG(5) 

CHARACTER ANSWER 

COMMON /BCR/ XTR, Y1TR, ADPATH, DATABLK, MYBUFSZ, NCHAN, AVG 
COMMON /CCSE/ TYME, DELTAFRQ, OUTFILE 
OUTFILE = 8 

* 

* OPEN FILES FOR OUTPUT OF DATA 

* 

OPEN (7JTLE=7usr/shafer/plate/data",STATUS="NEW") 

★ 

* CALCULATE HANNING WINDOW WEIGHTING FACTOR 

★ 

TPI = 8.0 * ATAN(l.O) 

PI = 4.0 * ATAN(l.O) 

TEMP = TPI / 1025.0 
SCL = SQRT(2.0/3.0) 

DO 401= 1,512 

WEGHT(I)= SCL * (1.0- COS (TEMP* FLO AT(I))) 

40 CONTINUE 

* 

* SET UP CONSTANTS FOR A/D CONVERTER SETUP AND BUFFERS 

* 

MYBUFSZ = 5000 
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NCHAN = 2 
FCHAN = 16 
SCHAN = 1 
GAIN = 0 
ADPATH = - 1 



* OPEN A/D CONVERTER PATH AND SET CONVERTER MODE 

* 

CALL MROPEN( ADPATH, "/dev/dacpO/adfl", RDONLY) 

CALL MRADENC(ADPATH, FCHAN, NCHAN, SCHAN, GAIN) 

CALL MRADMOD(ADPATH, 0, 0) 

* 

* OPEN CLOCK PATHD PULSES AND SET CLOCK MODE 

* 

* SET UP CLOCK FOR A/D CONVERTER 

* 

CLKPATH(l) = -1 
CLKPATH(2) = - 1 

CALL MROPEN (CLKPATH(l),7dev/dacp0/efclk0", RDWRT) 

CALL MROPEN (CLKPATH(2),7dev/dacpO/efclk4", RDWRT) 

FREQ = 5000. 

PRINT * ."SETTING UP CLOCKS" 

CALL MRCLK2 (CLKPATH(1),CLKPATH(2),0,0,FREQ,RFREQ,0,200000.,RWIDTH,2,0) 

CALL MRCLKTRIG (ADPATH, 2.CLKPATH) 

print Actual samp frequency per channel -- ",RFREQ 

NITEMS = 60000 

PRINT *, "NUMBER OF SAMPLES IS "NITEMS 

* 

* CALCULATE TIME RECORD LENGTH AND FREQUENCY RESOLUTION 

* 

RECLEN = 1024/RFREQ 

PRINT *,"TIME RECORD LENGTH IS ".RECLEN 

DFREQ= 1 /RECLEN 

DO 601= 1, 1024 

TYME(I) = I / RFREQ 

DELTAFRQ(I) = (I-1)*DFREQ 

60 CONTINUE 
* 

* GET NUMBER OF AVERAGES TO BE COMPUTED FROM USER 

* 

PRINT V INPUT NUMBER OF AVERAGES TO BE COMPUTED" 

READ * NOS AMP 

* 

* SET UP BUFFERS AND RELEASE THEM TO A/D PATH 

* 

DO 200 I = 1, 45001, MYBUFSZ 

CALL MRBUFID (ADPATH, DATABLK(I), MYBUFSZ*2) 

CALL MRBUFREL (ADPATH, DATABLK(I)) 

200 CONTINUE 
* 

* START DATA ACQUISITION 

* 

DO 250 L = 1 , NOSAMP 
65 PRINT *, L 
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PRINT V AWAITING TRIGGER ARM, TYPE ANY CHARACTER AND RETURN 
WHEN READY" 

READ *, ANSWER 
AVG(l) = 0 
AVG(2) = 0 

* 

* TAKE IN DATA RECORD 

* 

CALL MRXINQ (ADPATH, MYBUFSZ, NITEMS, ADSERVICE) MREATH.S, 30000 ) 
PRSTOPPED READING" 

* 

* DISPLAY DATA RECORD FOR CORRECTNESS 

* IF CORRECT CONTINUE, IF NOT DISCARD DATA SET AND TAKE ANOTHER 

* 

CALL PLOT1(1024, XTR,TYME,"DATA PREVffiW","VOLTAGE","TIME","pre.g") 

PRINT * , "Good data set?" 

READ *, ANSWER 
IF (ANSWER .EQ. "n") THEN 
GOTO 65 
END IF 

* 

* BEGIN DATA REDUCTION 

* 

DO 701= 1, 1024 
XX(I) = XTR(I) 

YY1(I) = Y1TR(I) 

70 CONTINUE 

* 

* WEIGHT TIME RECORD WITH HANNING WINDOW TO ENSURE PERIODICITY 

* REQUIREMENTS OF FFT ARE MET. XX, BEING AN IMPACT MEASUREMENT 

* IS PERIODIC BY DEFINITION AND REQUIRES NO WEIGHTING 

* 

DO 801= 1,512 
ITMP = 1025 - I 
YY1(I) = YY1(I) * WEGHT(I) 

YYl(ITMP) = YYIQTMP) * WEGHT(I) 

80 CONTINUE 
* 

* CALCULATE POWER SPECTRA INFORMATION 

* FOLLOWING PORTION TAKEN FROM CCSE PROGRAM 

* 

CALL FFT842(0, 1 024 ,XX , Y Y 1 .OUTFILE) 

GXX(l) = GXX(l) + 4.0 * XX(1)**2 
DO 90 K = 2, 513 
J = 1026 - K 

GXX(K) = GXX(K) + (XX(K)+XX(J))**2 + (YY1(K)-YY1(J))**2 
GYY(K) = GYY(K) + (YY1(K)+YY1(J))**2 + (XX(J)-XX(K))**2 
GXYRE(K) = GXYRE(K) + XX(K)*YY1(J) + XX(J)*YY1(K) 

GXYIM(K) =GXYIM(K) +XX(J)**2 +YY1(J)**2 -XX(K)**2 -YY1(K)**2 
90 CONTINUE 

GYY(l) = GYY(l) + 4.0* YY1(1)**2 
GXYRE(l) = GXYRE(l) + 2.0*(XX(1)*YY1(1)) 

GXYIM(l) = 0.0 
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CAUL ZERO(XTI, 1024) 

CALL ZERO(Y1TI,1024) 

* 

* CALCULATE FREQUENCY RESPONSE OF INPUT AND OUTPUT 

* SEPARATE OF THE POWER SPECTRA 

* 

DO 100 I = 1, 1024 
XX(I) = XTR(I) 

YY1(I) = Y1TR(I) 

100 CONTINUE 
* 

* WEIGHT TIME RECORD WITH HANNING WINDOW TO ENSURE PERIODICITY 

* REQUIREMENTS OF FFT ARE MET. XX, BEING AN IMPACT MEASUREMENT 

* IS PERIODIC BY DEFINITION AND REQUIRES NO WEIGHTING 

* 

DO 1101 = 1,512 
ITMP = 1025 - 1 
YY1(I) = YY1(I) * WEGHT(I) 

YYl(ITMP) = YY1(ITMP) * WEGHT(I) 

110 CONTINUE 

CALL FFT842(0,1024,XX,XTI,OUTFILE) 

CALL FFT842(0, 1024, YY 1 ,Y 1TI,0UTFILE) 

TEMP = 1/FLO AT(NOS AMP) 

DO 1201= 1, 513 
XFR(I) = XFR(I) + XX (I)* TEMP 
XFI(I) = XFI(I) + XTI(I)*TEMP 
Y1FR(I) = Y1FR(I) + YY1(I)*TEMP 
Y1FI(I) = Y1FI(I) + Y1TI(I)*TEMP 

120 CONTINUE 
* 

* TAKE NEXT DATA SET 

* 

250 CONTINUE 
* 

* WRITE DATA TO HARD DISK 

* 

DO 1301= 1,513 

WRITE(7,9990) GXX(I),GYY(I),GXYRE(I),GXYIM(I) 

130 CONTINUE 

9990 FORMAT(lX,G16.10,lx,G16.10,lx,G16.10,lxG16.10) 

* 

* WEIGHT POWER SPECTRA INFORMATION TO CORRECT FOR NUMBER OF DATA 

SETS 

* 

CALL CO HER( 1024 ,RFREQ,NOSAMP,GXX,GYY,GXYRE,GXYIM,l. 0,1.0) 

* 

* CALCULATE TRANSFER FUNCTION FROM POWER SPECTRA INFORMATION 

* 

DO 3001= 1,512 
K = 1025 - 1 

HFR(I) = GXYRE(I)/GXX(I) 

HFR(K) = HFR(I) 

HFI(I) = GXYIM(I)/GXX (I) 

HFI(K) = HFI(I) 

HF(I)= 10* ALOG(SQRT(HFR(I)**2+HFI(I)**2)) 
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HF(K)= HF(I) 

300 CONTINUE 
* 

* CALCULATE THE MOBILITY TRANSFER FUNCTION 

* 

DO 350 I = 2, 1024 

HMR(I) = HFR(I)/DELT AFRQ(I) 

HMI(I) = HFI(I)/DELTAFRQ(I) 

HM(I) = 10 * ALOG(SQRT(HMR(I)**2+HMI(I)**2)) 

350 CONTINUE 
* 

♦ CALCULATE THE FREQUENCY DOMAIN OF THE TIME SIGNALS 

* 

DO 400 I = 1, 513 

XFR(I) = SQRT(XFR(I)**2+XFT(I)**2) 

XFI(I) = 10* ALOG(XFR(I)) 

Y1FR(I) = SQRT(Y1FR(I)**2+Y1FI(I)**2) 

Y1FI(I) = 10 * ALOG(YlFR(I)) 

400 CONTINUE 
* 

* DISPLAY INFORMATION ON SCREEN AND PLOT IF DESIRED 

* 

CALL PL0T1(513, HER, DELTAFRQ,"H(0”, "Real", "Freq (Hz)","hfreq.g") 

CALL PLOT 1 (5 1 3.HFI.DELT AFRQ ,"Imag" ,"Freq (Hz)","pfreq.g") 

CALL PLOTl(513,HF,DELTAFRQ,"H(f)","db","Freq (Hz)”,"Hf.g”) 

CALL PLOTl(5 13, HMR,DELTAFRQ, "Mobility", "Real", "Freq (Hz)",”Hft.g") 

CALL PLOTl(5 13, HMI.DELT AFRQ, "Mobility", "Imag", "Freq (Hz)","Hft.g") 

CALL PLOTl(5 13, XFR.DELT AFRQ, ”X(f)", "Mag", "Freq (Hz)","xf.g") 

CALL PLOTl(513,XFI,DELTAFRQ,"X(f)","db","Freq (Hz)","logxf.g") 

CALL PLOT1(513,Y1FR,DELTAFRQ,"Y(0", "Mag”, "Freq (Hz)","yf.g") 

CALL PLOTl(5 13, Yin, DELT AFRQ, "Y(0","db", "Freq (Hz)","logyf.g") 

STOP 

END 

* 

* SUBROUTINE TO TAKE FILLED BUFFERS FROM A/D CONVERTER AND 

* CHECK TO SEE IF IMPACT HAS OCCURED. IF SO REMOVE DC COMPONENT 

* AND TRANSFER TO ARRAY FOR FURTHER REDUCTION AND RETURN BUFFER 

* TO READY BUFFER QUEUE AT A/D CONVERTER. IF NOT RETURN BUFFER 

* TO READY BUFFER QUEUE AT A/D CONVERTER. 

* 

SUBROUTINE ADSERVICE 



INTEGER 


K 


INTEGER 


INDEX 


INTEGER 


POWER 


INTEGER 


ADPATH 


INTEGER 


TEMP 


INTEGER 


NCHAN 


INTEGER 


MYBUFSZ 


INTEGER *2 


DATABLK( 100000) 


INTEGER 


GETINDEX 


PARAMETER 


(GETINDEX = 3) 


REAL 


DUMMY 


REAL 


AVG(5) 


REAL 


DUMMY 1 


REAL 


THOLD 
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REAL XTR(1024) 

REAL Y1TR(1024) 

COMMON /BCR/ XTR, Y1TR, ADPATH, DATABLK, MYBUFSZ, NCHAN, AVG 
TEMP = ADPATH 

CALL MRBUFGET (1, GETINDEX, INDEX) 

* 

* AVERAGE DATA SET TO DETERMINE DC OFFSET OF INFORMATION 

* 

IF (AVG(l) .EQ. 0.00) THEN 

DO 250 K = INDEX , INDEX + ( 99 * NCHAN), NCHAN 
AVG(l) = AVG(l) + FLOAT (DATABLK(K))/(2*20480) 

AVG(2) = AVG(2) + FLOAT(DATABLK(K+1))/(2*20480) 

250 CONTINUE 
END IF 

* 

* SET THRESHOLD TO 0.25V ABOVFSET OF HAMMER 

* 

THOLD = AVG(l) + 0.25 

DO 300 K = INDEX, INDEX + MYBUFSZ, NCHAN 
DUMMY =FLOAT(DATABLK(K))/(2*204.8) 

* 

* IF THRESHOLD IS EXCEEDED GO BACK 200 INDEXES AND RECORD 

* INFORMATION TO ARRAYS FOR FURTHER REDUCTION 

* 

IF ( DUMMY .GE. THOLD) THEN 
DO 350 I = -200 , 0, NCHAN 
J = K + I 
L = (1/2) + 101 

* 

* ENSURE BUFLAP IS CORRECT 

* 

IF (J .LE. 0) THEN 
J = J+ 100000 
END IF 

DUMMY = FLOAT(DATABLK(J)) 

DUMMY 1 = FLOAT(DATABLK(J+l)) 

XTR(L) = (DUMMY/(2*204 .8))- AVG( 1 ) 

Y1TR(L) = (DUMMYl/(2*204.8))-AVG(2) 

350 CONTINUE 
1 = 50 

* 

* READ IN REMAINING INFORMATION INTO DATA SETS FOR REDUCTION 

* 

DO 400 J = K,K + 975*NCHAN, NCHAN 
1 = 1+1 

DUMMY = FLOAT(DATABLK(J)) 

DUMMY 1 = FLOAT(DATABLK(J+l)) 

XTO(I) = (DUMMY/(2*204.8))-AVG(l) 

Y1TR(I) = (DUMMY 1/(2*204.8))- AVG(2) 

400 CONTINUE 
* 

* SET INDEX TO EXIT DATA TRANSFER LOOP 

* 

K = INDEX + MYBUFSZ + 1 
END IF 
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CONTINUE 



* RELEASE BUFFER TO READY BUFFER QUEUE FOR A/D CONVERTER 

* 

ADPATH = TEMP 

CALL MRBUFREL (TEMP, DAT AB LK(INDEX)) 

RETURN 
END 

SUBROUTINE COHER(NNN,ISR,NDSJP,GXX,GYY,GXYRE,GXYIM,SFX,SFY) 



xMAIN PROGRAM: A COHERENCE AND CROSS SPECTRAL ESTIMATION PROGRAxM 
AUTHORS: G. C. CARTER, J. F. FERRIE 

NAVAL UNDERWATER SYSTEMS CENTER 
NEW LONDON, CONNECTICUT 06320 
MODIFIED BY R.A. SHAFER, LT USN 

NAVAL POSTGRADUATE SCHOOL 
MONTERY CA 
APRIL 1987 

TO MEET THE REQUIREMENTS FOR ANALYSIS 
OF REAL TIME DATA 

INPUT: NNN IS THE NUMBER OF DATA POINTS PER SEGMENT 

4 < NNN < 1025 
ISR IS THE SAMPLING RATE 
NDSJP IS THE NUMBER OF DISJOINT SEGMENTS 
SFX IS THE SCALE FACTOR FOR THE INPUT DATA STORED IN 
THE XX ARRAY 

SFY IS THE SCALE FACTOR FOR THE INPUT DATA STORED IN 
THE YY ARRAY 



SPECIFICATION AND TYPE STATEMENTS 



REAL ISR 

REAL XX(1024), YY(1024) 

REAL GXX(513), GYY(513), GXYRE(513), GXYIM(513) 
REAL WEGHT(513), PHI(513) 

REAL LINE(50) 

REAL DELTAFRQ(1024),TYME(1024) 

EQUIVALENCE (WEGI(l)) 



COMMON /CCSE/ TYME, DELTAFRQ, OUTFILE 



NNN IS THE NUMBER OF DATA POINTS PER SEGMENT 

ISR IS THE SAMPLING RATE 

NDSJP IS THE NUMBER OF DISJOINT SEGMENTS 

SFX AND SFY ARE SCALE FACTORS FOR THE INPUT DATA 



NFFTS = NDSJP 
SMALL = .lE-20 



CALCULATE CONSTANTS 



TPI = 8.0*ATAN(1.0) 
DEG = 360.0/TPI 
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VARX = 0.0 
VARY = 0.0 
DT = 1.0/ FLOAT(ISR) 

SF = SQRT(ABS(SFX*SFY)) 

NPFFT = NNN 

NNNP1 = NNN + 1 

NNND2 = NNN/2 

NND21 = NNND2 + 1 

NP2 = NPFFT + 2 

ND2 = NPFFT/2 

ND2P1 = ND2 + 1 

DF = 1 .0/(DT* FLOAT (NPFFT)) 

FNYQ = FLOAT(ISR)/2.0 
CONST = 0.25*DT/FLOAT(NNN) 

FLOW = 0.0 
FfflGH = FNYQ 
ISTRT = DFIX(FLOW/DF) + 1 
ISTOP = DFIX(FfflGH/DF) + 1 
C 

C NORMALIZE ESTIMATES 
C 

FNSG = FLOAT(NFFTS) 

OFNSG = 1.0/FNSG 
TEMPI = CONST*OFNSG*SFX 
TEMP2 = CONST*OFNSG*SFY 
TEMP4 = CONST*OFNSG*SF 
TEMP3 = 2.0*TEMP4 
DO 90 K=UND2P1 
GXX(K) = GXX(K)*TEMP1 
GYY(K) = GYY(K)*TEMP2 
GXYRE(K) = GXYRE(K)*TEMP3 
GXYIM(K) = GXYIM(K)*TEMP4 
90 CONTINUE 
VARX = 0.0 
VARY = 0.0 
DO 100 K=1,ND2P1 
VARX = VARX + GXX(K) 

VARY = VARY + GYY(K) 

100 CONTINUE 

VARX = VARX*DF*2.0/SFX 
VARY = VARY*DF*2.0/SFY 
C 

C CONVERT GXX TO DB AND PLOT 
C 

DO 110 1=1 J4D2P1 
XX(I) = GXX (I) 

PHI(I)= 10.0*ALOG10(AMAX1(GXX(I),SMALL)) 

110 CONTINUE 

CALL PLOTl(ND2Pl r PHI,DELTAFRQ,"Gxx (f)","db","FREQ (Hz)", ,, loggxx.g") 
C 

C CONVERT GYY TO DB AND PLOT 
C 

DO 140 I=1,ND2P1 
XX(I) = GYY(I) 

PHI (I) = 1 0.0* ALOG 1 0( AMAX 1 (GYY (I) ,S MALL)) 
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140 CONTINUE 

CALL PL0T1(ND2P1, PHI, DELTAFRQ, "Gyy(0","db","FREQ (Hz)","gyy.g") 
COMPUTE CROSS SPECTRUM AND MAGNITUDE SQUARED COHERENCE 
DO 230 K=1ND2P1 

PHI(K) = GXYR£(K)**2 + GXYIM(K)**2 
XX(K) = SQRT(ABS(PHI(K)/(GXX(K)*G YY (K)))) 

30 CONTINUE 

CALL PL0T1(ND2P1, XX, DELTAFRQ, "Coherence", " ", "Freq (Hz)","coher.g") 

TERMINATE PROGRAM 

RETURN 

END 



SUBROUTINE: ZERO 

THIS SUBROUTINE STORES IN A FLOATING POINT ARRAY 



SUBROUTINE ZERO(ARRAY, NUMBR) 

INPUT: ARRAY = AN ARRAY OF FLOATING POINT VALUES TO BE 
ZERO FILLED 

NUMBR = NUMBER OF ARRAY VALUES 

DIMENSION ARRAY(l) 

C 

DO 10 K=l, NUMBR 
ARRAY(K) = 0.0 
10 CONTINUE 
RETURN 
END 
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APPENDIX D 



STEADY STATE VIBRATION DATA ACQUISITION AND REDUCTION 

PROGRAM 

This Appendix list the computer program utilized to acquire the data for steady state 
vibrations analysis. All routines were written by the author except where noted. For 
further explanation of the data acquisition specific routines the reader is encouraged to refer 
to the MASSCOMP documentation "Data Acquisition User's Manual". Figure 48 is the 
program flow diagram. Sections of the program follow the general outline presented in the 
flow diagram. Specific subroutines for plotting the information are not shown in this 
Appendix. 



♦ THIS PROGRAM IS THE MAIN PROGRAM FOR DATA ACQUISITION 

* FOR STEADY STATE ANALYSIS OF A VIBRATING BEAM 



* 07 MAY 1987 

* 



INTEGER 
INTEGER 
INTEGER 
INTEGER 
REAL 
REAL 
INTEGER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
INTEGER *2 
INTEGER 
EQUIVALENCE 
EXTERNAL 
INTEGER 
INTEGER 
INTEGER 
INTEGER 



MYBUFSZ, N BUFFS 
NEARFREQ, RDONLY 
RDWRT, SQUARE 
STARTLO, TIMEOUT 
FREQ, WIDTH 
PI 

OUTFDLE 

( NBUFFS = 5 ) 

( NEARFREQ = 0 ) 

( RDONLY = 1 ) 

( RDWRT = 0 ) 

( SQUARE = 4 ) 

( STARTLO = 0 ) 

( TIMEOUT = 30000 ) 

( WIDTH = 0.0 ) 
DATABLK(40960) 
DUMMY 

(DATABLK, DUMMY) 
ADSERVICE 
ADPATH 
CLKPATH(2) 
STATWDS(2) 

I 
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INTEGER NCHAN 

INTEGER FCHAN 

INTEGER SCHAN 

INTEGER GAIN 

INTEGER INDEX 

INTEGER POWER 

REAL TPI 

REAL WEGHT(1025) 

REAL TYME(2048) 

REAL RECLEN 

REAL RFREQ 

REAL RWIDTH 

REAL XTR(2048), Y1TR(2048) 

REAL Y1TI(2048) 

REAL XX(2048) 

REAL YY 1(2048) 

REAL XFR(2048),XFI(2048) 

REAL Y1FR(2048),Y1FT(2048) 

REAL GXX(1025) 

REAL GYY(1025) 

REAL GXYRE(1025) 

REAL GXYIM(1025) 

REAL HWEN(2048) 

REAL DELTAFRQ(2048) 

REAL DFREQ 

REAL MAG1(1025) 

REAL MAG2(1025) 

REAL MAG3(1025) 

REAL HFR(2048),HFI(2048) 

REAL HR(2048),HI(2048) 

REAL AVG(5) 

REAL PI2, TEMP 

INTEGER AVGSAMP 

CHARACTER ANSWER 

COMMON /BCR/ XTR, Y1TR, ADPATH, DATABLK, MYBUFSZ, NCHAN, AVG 
COMMON /CCSE/ TYME, DELTAFRQ, OUTFILE 

* 

* OPEN DATA FILES FOR DATA STORAGE 

* 

OUTFILE = 8 

OPEN (7/ile="/usr/shafer/beam/simuldata") 

* 

* CALCULATE WINDOWS AND CONSTANTS 

* 

TPI = 8.0 * ATAN(l.O) 

PI = 4.0 * ATAN(l.O) 

TEMP = TPI / 2049.0 
SCL = SQRT(2.0/3.0) 

DO 401 = 1,1024 

WEGHT(I)= SCL * (1.0- COS(TEMP*FLOAT(I))) 

40 CONTINUE 

PI2 = 2.0 * ATAN(l.O) 

TEMP = PI2/ 512.0 
DO 46 I = 2, 1024 
HWIN(I) = 2.0 
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46 CONTINUE 

DO 47 I = 1025,2048 
HWEN(I) = 0.0 

47 CONTINUE 
HWTN(l) = 1.0 

* 

* SET UP A/D CONVERTER AND BUFFERS FOR TEMPORARY STORAGE OF 

DATA 

4c 

MYBUFSZ = 4096 
NCHAN = 2 
FCHAN = 16 
SCHAN = 1 
GAIN = 0 
AD PATH = -1 

CALL MROPEN(ADPATH, "/dev/dacpO/adfl", RDONLY) 

CALL MRADINC(ADPATH, FCHAN, NCHAN, SCHAN, GAIN) 

CALL MRADMOD(ADPATH, 0, 0) 

* INDENTIFY BUFFERS AND RELEASE THEM TO READY BUFFER QUEUE 

* 

DO 200 I = 1, 36865, MYBUFSZ 
CALL MRBUFID (ADPATH, DATABLK(I), MYBUFSZ*2) 

CALL MRBUFREL (ADPATH, DATABLK(I)) 

200 CONTINUE 

4c 

* SET UP CLOCK FOR A/D CONVERTER 

4c 

CLKPATH(l) = -1 
CLKPATH(2) = -1 

CALL MROPEN (CLKPATH(l),"/dev/dacpO/efclkO", RDWRT) 

CALL MROPEN (CLKPATH(2),"/dev/dacpO/efclk4", RDWRT) 

FREQ = 1000. 

PRINT *, "SETTING UP CLOCKS" 

CALL 

MRCLK2(CLKPATH(1), CLKPATH(2),0,0,FREQ,RFREQ,0, 200000., RWIDTH, 2,0) 

CALL MRCLKTRIG (ADPATH, 2.CLKPATH) 

4c 

* COMPUTE FREQUENCY RESOLUTION AND FILL FREQUENCY ARRAY 

* DISPLAY RESOLUTION ON SCREEN 

* 

RFREQ = RFREQ 

PRINT *, "Actual sampling frequency per channel - ", RFREQ 
PRINT V INPUT NUMBER OF AVERAGES TO BE COMPUTED" 

READ *,AVGSAMP 
NITEMS = 4096 
RECLEN = 2048/RFREQ 

PRINT VTIME RECORD LENGTH IS PER CHANNEL IS "RECLEN 
DFREQ= 1 /RECLEN 

PRINT FREQUENCY RESOLUTION IS ",DFREQ 

DO 601 = 1,2048 
TYME(I) = I / RFREQ 
DELTAFRQ(I) = (I-1)*DFREQ 
60 CONTINUE 
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* MAIN LOOP FOR DATA ACQUSITION AND PRELIMINARY REDUCTION 

* 

DO 250 L = 1 , AVGSAMP 

CALL MRXINQ (ADPATH, MYBUFSZ, NITEMS, ADSERVICE) 

CALL MREVWT (ADPATH, STATWDS, 0 ) 

PRINT *, AVG(3) 

DO 701= 1,2048 
XX(I) = XTR(I) 

YY1(I) = Y1TR(I) 

70 CONTINUE 
* 

* WEIGHT TIME RECORDS WITH HANNING WINDOW 

* 

DO 80 I = 1, 1024 
ITMP = 2049 - 1 
XX(I) = XX(I) * WEGHT(I) 

YY1(I) = YY1(I) * WEGHT(I) 

XX(ITMP) = XX(ITMP) * WEGHT(I) 

YY 1 (ITMP) = YY1(ITMP) * WEGHT(I) 

80 CONTINUE 
* 

* COMPUTE FORWARD FFT 

* 

CALL FFT842(0, 2048, XX, YYl.OUTFILE) 

* 

* COMPUTE SPECTRA 

* 

GXX(l) = GXX(l) + 4.0*XX(1)**2 
DO 90 K=2, 1025 
J = 2050 - K 

GXX(K) = GXX(K) + (XX(K)+XX(J))**2 + (YY1(K)-YY1(J))**2 
GYY(K) = GYY(K) + (YY1(K)+YY1(J))**2 + (XX(J)-XX(K))**2 
GXYRE(K) = GXYRE(K) + XX(K)*YY1(J) + XX(J)*YY1(K) 

GXYIM(K) = GXYIM(K) +XX(J)**2 +YY1(J)**2 -XX(K)**2 -YY1(K)**2 
90 CONTINUE 

GYY(l) = GYY(l) + 4.0*YY1(1)**2 
GXYRE(l) = GXYRE(l) + 2.0*(XX(1)*YY1(1)) 

GXYIM(l) = 0.0 



250 CONTINUE 
* 

* DATA COLLECTION COMPLETE, COMPLETE REDUCTION AND DISPLAY 

RESULTS 

* 

CALL COHER1 (2048 JRPREQAVGS AMP, GXX,GYY,GXYRE,GXYIM, 1.0, 1.0) 

* 

* CALCULATE FREQUENCY RESPONSE (MOBILITY) 

* 

DO 600 I = 2, 1024 
K = 2049 - 1 

HFR(I) = (GXYRE(I)/GXX(I))/DELTAFRQ(I) 

HFR(K) = HFR(I) 



79 



uuuuuuuuuuuuuuuuuuuuuuuu 



HFI(I) = (GXYIM(I)/GXX(I))/DELTAFRQ(I) 

HFI(K) = HFI(I) 

600 CONTINUE 

HFR(l) = GXYRE(1)/GXX(1) 

HFR(2048) = HFR(l) 

HFI(l) = GXYIM(1)/GXX(1) 

HFI(2048) = HFI(l) 

* 

* RECORD POWER SPECTRA INFORMATION TO DISK 

* 

DO 700 1 = 1, 1025 

WRITE (7,9990) GXX(I),GYY(I),GXYRE(I),GXYIM(I) 

700 CONTINUE 

9990 FORMAT (IX, G16.8,G16.8,G16.8,G16.8) 

* 

* PLOT RESULTS 

* 

CALL PLOT1 (1024 ,HFR,DELTAFRQ,"H(f)", "Real", "Freq (Hz)Y'Hfr.g") 

CALL PLOT1 (1024, HFI,DELTAFRQ,"H(f)","Imag", "Freq (H 2 )","Hfi.g") 

STOP 

END 

SUBROUTINE COHERl(NNN,ISR,NDSJP,GXX,GYY,GXYRE,GXYIM,SFX,SFY) 



MAIN PROGRAM: A COHERENCE AND CROSS SPECTRAL ESTIMATION PROGRAM 
AUTHORS: G. C. CARTER, J. F. FERRIE 

NAVAL UNDERWATER SYSTEMS CENTER 
NEW LONDON, CONNECTICUT 06320 
MODIFIED BY R.A. SHAFER, LT USN 

NAVAL POSTGRADUATE SCHOOL 
MONTERY CA 
APRIL 1987 

TO MEET THE REQUIREMENTS FOR ANALYSIS 
OF REAL TIME DATA 

INPUT: NNN IS THE NUMBER OF DATA POINTS PER SEGMENT 

4 < NNN < 1025 
ISR IS THE SAMPLING RATE 
NDSJP IS THE NUMBER OF DISJOINT SEGMENTS 
SFX IS THE SCALE FACTOR FOR THE INPUT DATA STORED IN 
THE XX ARRAY 

SFY IS THE SCALE FACTOR FOR THE INPUT DATA STORED IN 
THE YY ARRAY 



SPECIFICATION AND TYPE STATEMENTS 
REAL ISR 

REAL XX(2048), YY(2048) 

REAL GXX(1025), GYY(1025), GXYRE(1025), GXYIM(1025) 
REAL WEGHT(1025), PHI(1025) 

REAL LINE(50) 

REAL DELT AFRQ(204 8) ,TY ME(204 8) 

EQUIVALENCE (WEGHT(l)FHI(l)) 

COMMON /CCSE/ TYME, DELTAFRQ, OUTFILE 
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c 

C NNN IS THE NUMBER OF DATA POINTS PER SEGMENT 
C ISR IS THE SAMPLING RATE 
C NDSJP IS THE NUMBER OF DISJOINT SEGMENTS 
C SFX AND SFY ARE SCALE FACTORS FOR THE INPUT DATA 

NFFTS = NDSJP 
SMALL = .lE-20 



CALCULATE CONSTANTS 

TPI = 8.0*ATAN(1.0) 

DEG = 360.0/TPI 
VARX = 0.0 
VARY = 0.0 
DT= 1.0/FLOAT(ISR) 

SF = SQRT(ABS(SFX*SFY)) 

NPFFT = NNN 

NNNP1 = NNN + 1 

NNND2 = NNN/2 

NND21 = NNND2 + 1 

NP2 = NPFFT + 2 

ND2 = NPFFT/2 

ND2P1 = ND2 + 1 

DF = 1 .0/(DT* FLOAT (NPFFT)) 

FNYQ = FLOAT(ISR)/2.0 

CONST = 0.25*DT/FLOAT(NNN) 

FLOW = 0.0 

FHIGH = FNYQ 

ISTRT = IFIX(FLOW/DF) + 1 

ISTOP = ffIX(FHIGH/DF) + 1 

NORMALIZE ESTIMATES 

FNSG = FLOAT(NFFTS) 

OFNSG = 1.0/FNSG 
TEMPI = CONST*OFNSG*SFX 
TEMP2 = CONST*OFNSG*SFY 
TEMP4 = CONST*OFNSG*SF 
TEMP3 = 2.0*TEMP4 
DO 90 K=1 JND2P1 
GXX(K) = GXX(K)*TEMP1 
GYY(K) = GYY(K)*TEMP2 
GXYRE(K) = GXYRE(K)*TEMP3 
GXYIM(K) = GXYIM(K)*TEMP4 
90 CONTINUE 
VARX = 0.0 
VARY = 0.0 
DO 100 K=1,ND2P1 
VARX = VARX + GXX(K) 

VARY = VARY + GYY(K) 

100 CONTINUE 

VARX = VARX*DF*2.0/SFX 
VARY = VARY*DF*2.0/SFY 
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C CONVERT GXX TO DB AND PLOT 
C 

DO 110 I=1,ND2P1 
XX(I) = GXX(I) 

PHI (I) = 10.0* ALOG10(AMAX1(GXX(I), SMALL)) 

110 CONTINUE 

CALL PL0T1(ND2P1,PHI, DELTAFRQ, "Gxx (0","db";’FREQ (Hz)","loggxx.g") 

C 

C CONVERT GYY TO DB AND PLOT 
C 

DO 140 I=1,ND2P1 
XX(I) = GYY(I) 

PHI(I)= 1 0.0* ALOG10(AMAX1(GYY(I), SMALL)) 

140 CONTINUE 

CALL PLOTl(ND2Pl, PHI, DELTAFRQ, "Gyy(0' , , , ’db","FREQ (Hz)”,"gyy.g") 

C 

C COMPUTE AND DISPLAY PHASE FROM AVERAGED GXYRE AND GXYIM SPECTRUM 
C 

GXYIM(l) = 0.0 
PHI(l) = 0.0 
DO 200 K=2,ND2P1 
XXK = GXYRE(K) 

IF (XXK) 190, 170, 190 
170 IF (GXYIM(K)) 190, 180, 190 
180 XXK =10 

190 PHI(K) = DEG*ATAN2(GXYIM(K),XXK) 

200 CONTINUE 
C 

C PLOT PHASE FROM -PHLIM TO PHLIM 
C 

PHLIM = 1800.0 
DO 210 K=2,ND2P1 
X = PHI(K) - PHI(K-l) 

PHI(K) = PH3(K) - SIGN(360.,X)*AINT(0.5+ABS(X)/360.0) 

IF (PHI(K).GT.PHLIM) PHI(K) = PHI(K) - PHLIM 
IF (PHI(K).LT.(-PHLIM)) PHI(K) = PHI(K) + PHLIM 
210 CONTINUE 
C 

C COMPUTE CROSS SPECTRUM AND MAGNITUDE SQUARED COHERENCE 
C 

DO 230 K=1,ND2P1 

PHI(K) = GXYRE(K)**2 + GXYIM(K)**2 
XX(K) = SQRT(ABS(PHI(K)/(GXX(K)*GYY(K)))) 

230 CONTINUE 

CALL PLOTl(ND2Pl,XX, DELTAFRQ, "Coherence", " ”, "Frcq (hz)","coher.g") 

C 

C CONVERT GXY TO DB AND PLOT 
C 

DO 250 I=1,ND2P1 

PHI(I) = 5.0* ALOG10( AMAX1(PHI(I). SMALL)) 

250 CONTINUE 

C 

C TERMINATE PROGRAM 
C 

RETURN 
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END 



SUBROUTINE: ZERO 

THIS SUBROUTINE STORES ZEROES IN A FLOATING POINT ARRAY 



SUBROUTINE ZERO(ARRAY, NUMBR) 

INPUT: ARRAY = AN ARRAY OF FLOATING POINT VALUES TO BE 
ZERO FILLED 

NUMBR = NUMBER OF ARRAY VALUES 

DIMENSION ARRAY(l) 

C 

DO 10 K=1 .NUMBR 
ARRAY(K) = 0.0 
10 CONTINUE 
C 

RETURN 

END 

* 

* SUBROUTINE FOR DATA ARRICE, TRANSFER INTO WORKING ARRAY 

* 

SUBROUTINE ADSERVICE 
INTEGER K 

INTEGER INDEX 

INTEGER POWER 

INTEGER ADPATH 

INTEGER TEMP 

INTEGER NCHAN 
INTEGER MYBUFSZ 

INTEGERS DATABLK(40960) 

INTEGER GETINDEX 

PARAMETER (GETINDEX = 3) 

REAL DUMMY 

REAL AVG(5) 

REAL DUMMY 1 

REAL THOLD 

REAL XTR(2048) 

REAL Y1TR(2048) 

COMMON /BCR/ XTR, Y1TR, ADPATH, DATABLK, MYBUFSZ, NCHAN, AVG 

TEMP = ADPATH 

AVG(3) = AVG(3) + 1 

AVG(l) = 0 

AVG(2) = 0 

★ 

* GET DONE BUFFER FROM DONE BUFFER QUEUE OF A/D CONVERTER 

* 

CALL MRBUFGET (1, GETINDEX, INDEX) 

★ 

* AVERAGE ENTIRE ARRAY 

* 

DO 250 I = INDEX, INDEX + MYBUFSZ , 2 

AVG(l) = AVG(l) + FLOAT(DATABLK(I))/(2.* 419430.4) 

AVG(2) = AVG(2) + FLOAT(DATABLK(I+l))/(2. * 419430.4) 
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250 CONTINUE 
L = 0 

* 

* READ IN ARRAY INTO WORKING ARRAY, CONVERTING 12 BIT RESOLUTION 

* INTO CORRECT FORMAT (2*204.8 IS CORRECTION FACTOR) 

* 

DO 300 K = INDEX, INDEX + MYBUFSZ, 2 
L = L + 1 

DUMMY = FLOAT(DATABLK(K)) 

DUMMY 1 = FLOAT(DATABLK(K+l)) 

XTR(L) = (DUMMY/(2.*204.8))-AVG(l) 

Y1TR(L) = (DUMMY 1/(2. *204.8))-AVG(2) 

300 CONTINUE 
* 

* RELEASE BUFFER TO READY BUFFER QUEUE OF A/D CONVERTER 

* 

ADPATH = TEMP 

CALL MRBUFREL (TEMP, DATABLK(INDEX)) 

RETURN 

END 
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Figure 48: Steady State Vibrations Data Acquisition and Reduction Program Flow 

Diagram 
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